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Abstract  -  In  this  paper  we  present  a  recursive  Bayesian 
solution  to  the  problem  of  joint  tracking  and  classification  of  a 
target  modeled  at  a  distance  by  an  equivalent  magnetic  dipole. 
Tracking/classification  of  a  magnetic  dipole  from  noisy 
magnetic  field  measurements  involves  the  modeling  of  a 
non-linear  non-Gaussian  system.  This  system  allows  for 
complications  due  to  multiple  directions  of  arrival  and  target 
maneuver.  The  determination  of  target  position,  velocity  and 
magnetic  moment  is  formulated  as  an  optimal  stochastic 
estimation  problem,  which  could  be  solved  using  a  sequential 
Monte  Carlo  based  approach  known  as  the  particle  filter.  In 
addition  to  the  conventional  particle  filter,  the  proposed 
tracking  and  classification  algorithm  uses  the  unscented 
Kalman  filter  (UKF)  to  generate  the  transition  prior  as  the 
proposal  distribution. 

I.  INTRODUCTION 

Precise  determination  of  target  motion  parameters,  i.e. 
position,  velocity,  and  target  classification,  are  primary 
concerns  in  automated  surveillance  systems.  This  paper 
presents  the  use  of  a  sequential  Monte  Carlo  based 
statistical  signal  processing  method,  known  as  particle  filter, 
for  tracking  and  classifying  a  magnetic  target. 

A  target  containing  ferromagnetic  material  can  be 
adequately  modeled  at  a  distance  by  an  equivalent  magnetic 
dipole  moment.  This  magnetic  target  can  be  observed  by 
means  of  tri-axial  magnetometers  that  measure  the  variation 
of  the  magnetic  field  components  as  a  function  of  time  as  it 
passes.  Of  interest  is  the  inverse  problem  of  the 
determination  of  the  position  and  magnetic  parameters  of 
the  target  at  time  step  k  from  its  magnetic  signature 
collected  up  to  and  including  time  k. 

One  approach  to  solve  this  problem  makes  use  of  the 
recursive  Bayesian  estimation  (filtering)  technique.  The 
problem  is  formulated  in  state-space  form  where  the  state 
variables  are  the  position,  velocity  and  magnetic  moment  of 
the  target.  Let  define  xk  as  the  state  of  the  system  at  time 
step  k,  and  Zj:k  =  {zj,  z2,  ...,  zk}  as  the  observation 
(measurements)  history  of  a  system  from  time  1  to  k. 
Because  of  either  noise  in  the  state  evolution  process  or 
uncertainty  as  to  the  exact  nature  of  the  process  itself,  the 
state  vector  xk  is  generally  regarded  as  a  random  variable.  In 
the  Bayesian  filtering  technique,  one  attempts  to  construct 
an  estimate  of  the  posterior  probability  density  function 
(pdf),  p(xk  |  Zj:k).  Since  all  information  provided  by  zkkis 
conveyed  by  the  posterior  density,  it  may  be  said  to  be  the 
complete  solution  to  the  estimation  problem.  A  recursive 
Bayesian  algorithm  imposes  the  constrain  that  the  estimate 
of  p(xk  |  z1:k)  should  be  generated  solely  from  the  previous 


posterior  density,  p(xk_i  |  Zi:k_j),  and  the  most  recent 
measurement  zk.  In  this  way,  it  is  not  necessary  to  store  the 
complete  data  set  or  to  reprocess  existing  data  when  a  new 
measurement  becomes  available. 

The  recursive  propagation  of  the  posterior  density  is 
only  a  conceptual  solution  that  can  be  determined 
analytically  only  in  a  restrictive  set  of  cases.  When  the 
analytical  solution  is  intractable,  a  Monte  Carlo  based 
approach  to  recursive  Bayesian  filtering  called  the  particle 
filter,  is  one  method  that  approximates  the  optimal  Bayesian 
solution.  In  the  Monte  Carlo  method,  a  set  of  random 
samples  (particles)  are  drawn  from  a  target  distribution  such 
as  p(x  |  z).  In  general,  this  distribution  is  not  known.  We  will 
use  q(xk  |  z1:k)  *  p(xk  |  z1:k)  to  denote  a  proposal  distribution 
from  which  samples  can  be  drawn.  The  main  drawback  of 
the  conventional  particle  filter  is  that  it  uses  transition  prior, 
p(xk  |  xk_i),  as  the  proposal  distribution.  The  transition  prior 
does  not  take  into  account  current  observation  data.  To 
overcome  this  difficulty,  the  unscented  Kalman  filter  (UKF) 
was  proposed  to  generate  better  proposal  distributions  by 
taking  into  consideration  the  most  recent  observation. 


II.  UNSCENTED  PARTICLE  FILTER 

A.  Recursive  Bayesian  Estimation 

The  tracking  problem  requires  estimation  of  the  state 
vector  (target  co-ordinates,  velocity,  magnetic  moment)  of  a 
system  that  changes  over  time  using  a  sequence  of  noisy 
measurements  (observations)  made  on  the  system.  For  the 
specific  application  regarding  this  study,  the  target 
dynamics  (the  system  model)  is  described  by  a  linear 
equation,  f(«),  while  the  system  observation  (the 
measurement  model)  equation,  h(«),  is  highly  non-linear. 
We  assume  that  these  models  are  available  in  a  probabilistic 
form: 


P(xk  1  Vi) : 

xk  =f(xH-Vl) 

(1) 

P(2t  lAt): 
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where  xk  is  the  nx-dimensional  state  vector  of  the  system  at 
time  step  k,  zk  is  the  nz-dimensional  observation  vector,  and 
vk  and  wk  are  vectors  representing  the  process  and 
measurement  noise,  respectively.  They  have  the  dimensions 
nv  and  nw.  It  is  assumed  that  the  noise  vectors  are  i.i.d.  and 
independent  of  current  and  past  states. 

From  the  Bayesian  perspective,  it  is  required  to  estimate 
p(xk  |  z1:k)  assuming  that  the  pdf  at  time  (£-1),  p(xk_!  |  z,  :k-i), 
is  available.  The  first  step  in  this  process  is  called  prediction 
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and  makes  use  of  equation  (1),  which  is  assumed  to  describe 
a  Markov  process  of  order  one: 


P(*k  I  Zl:*-1 )  =  f  P(*k  I  XA-1 )  P(*k- 1  I  Z1:A- 1 )  ^XA-1  (3) 


The  second  step,  the  measurement  update,  uses  the  most 
recent  observation  to  produce  the  desired  pdf  via  Bayes’ 
rule: 


P(*k\zu)  = 


p(zk\xk)  p(xk\zvk_,) 


P(Zk  lZbt-l) 

P( Zk  I  Z1:A-1  )=\p(Zk  I  XA- )  P(Xk  I  Zl:A- 1 )  rfXA 


(4) 


where  the  second  equation  is  the  normalization  constant. 
Once  the  posterior  pdf  is  determined,  it  is  straightforward 
conceptually  to  produce  any  desired  statistic  of  xk.  For 
instance,  the  minimum  mean-square  error  (MMSE) 
estimate  of  the  current  state  could  be  found  by  computing 
the  conditional  mean: 


x,  = 


=  {XA  P(*k 


. )  r/x  /. 


(5) 


The  conditional  covariance  matrix  is  obtained  in  a 
similar  way: 

Pi  =  {  (XA  -  4  )  (X*  -  4  f  P^k  I  Z  Ik  )  dxk  (6) 


unknown  distribution,  p,  by  a  set  of  properly  weighted 
particles  drawn  from  a  known  distribution,  q.  In  this  way, 
the  difficult  problem  of  distribution  estimation  is  converted 
to  an  easy  problem  of  weight  estimation.  The  exact  form  of 
the  proposal  distribution  q  is  a  critical  issue  in  designing  the 
particle  filter  and  is  usually  approximated  to  facilitate  easy 
sampling. 

A  numerical  approximation  to  the  recursive  Bayesian 
filtering  method  given  by  the  equation  (3)  and  (4)  is  the 
following  algorithm: 


1.  Initialization:  sample  N  particles  xk(l),  i  =  1,  2,  ...,  N, 
from  the  proposal  distribution.  The  proposal 
distribution  can  be  the  transition  prior  as  used  in  the 
conventional  particle  filters,  or  more  advanced 
distributions  like  the  one  used  in  this  study. 

2.  Measurement  update:  update  the  importance  weights. 
The  Bayesian  sequential  importance  sampling  (SIS) 
procedure  gives  a  recursive  calculation  of  the 
normalized  weight  [1]: 
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In  general,  the  recursive  propagation  of  the  posterior 
density  cannot  be  determined  analytically  because  the 
integrals  in  (3)  and  (4)  do  not  have  closed-form  solutions. 
Solutions  do  exist  in  a  restrictive  set  of  cases.  For  example, 
if  f(»)  and  h(»)  are  linear  functions  and  if  Gaussian 
distributions  are  assumed  for  x,  v,  and  w,  the  estimation  of 
states  is  reduced  to  the  well-known  Kalman  filter. 

The  problem  of  tracking  a  magnetic  dipole  does  not 
satisfy  the  original  Kalman  filter  requirements  because  the 
system  observation  is  non-linear.  Moreover,  because  the 
target  can  approach  the  sensors  from  any  direction  and  can 
maneuver  at  any  time,  the  true  posterior  density  is 
multi-modal  and  a  Gaussian  description  will  be  inaccurate. 

In  order  to  deal  with  non-linear  systems  and/or 
non-Gaussian  reality,  two  categories  of  techniques  have 
been  developed:  parametric  and  non-parametric.  The 
parametric  techniques  are  based  on  improvements  of  the 
Kalman  filter.  These  filters  (for  example,  extended  and 
unscented  Kalman  filters)  can  handle  non-linear  equations, 
but  they  implicitly  approximate  the  posterior  density  as 
Gaussian.  The  non-parametric  techniques  are  based  on 
Monte  Carlo  simulations  and  are  the  subject  of  the  present 
study.  These  filters  assume  no  functional  form,  but  instead 
use  a  set  of  random  samples  (particles)  to  estimate  the 
posteriors.  The  advantage  is  that  the  particle  filters  can 
accommodate  simultaneous  alternative  hypotheses  that  can 
describe  a  multi-modal  distribution  well. 

B.  Particle  Filter  Implementation 

The  basic  idea  of  the  Monte  Carlo  based  approach  to  an 
intractable  Bayesian  filtering  case  is  to  approximate  an 


As  an  approximation  to  (5)  take: 


i=l,N 

3.  Re-sampling  is  a  necessary  step  introduced  in  particle 
filtering  algorithms  to  reduce  the  degeneration  of 
samples.  In  practice  it  was  noticed  that,  after  a  few 
iterations,  one  of  the  importance  weights  tends  to  one, 
while  the  others  become  zero.  To  avoid  the  degeneracy, 
the  sampling  importance  re-sampling  (SIR)  method 
selects  N  samples  with  replacement  from  the  set  xk(l), 
where  the  probability  to  take  sample  ‘i’  is  wk(1).  Then  set 
wk(i)  =  l/N,  i  =  1,  2,  . . .,  N. 

4.  Prediction:  assuming  that  the  probability  of  the 
process  noise  is  known,  use  equation  (1)  to  simulate 
xk+1(i),i=l,2,  ...,  N. 

5.  Set  k  =  k +1,  and  iterate  to  item  2. 

C.  The  Unscented  Kalman  Filter 

As  mentioned,  the  deficiency  of  the  sequential 
importance  sampling  (SIS)  approximation  is  that  the 
proposal  distribution  may  be  very  different  from  the 
posterior  distribution,  especially  if  using  the  transition  prior 
as  the  proposal  distribution.  An  improved  proposal 
distribution  must  incorporate  the  current  observation  data 
with  the  optimal  Gaussian  approximation  of  the  state. 

In  a  previous  study  [2]  on  the  magnetic  dipole  tracking 
application,  it  was  shown  that  the  unscented  Kalman  filter 


(UKF)  is  the  best  Kalman  filter  for  the  non-linear  systems. 
The  UKF  is  so  named  because  it  implements  the  Kalman 
recursion  using  the  sample  points  provided  by  the  unscented 
transform.  The  unscented  transform  deterministically 
generates  a  set  of  points  that  have  a  certain  mean  and  sample 
covariance.  The  non-linear  function  is  then  applied  to  each 
of  the  sample  points,  yielding  a  transformed  sample  from 
which  the  predicted  mean  and  covariance  are  calculated. 
The  estimate  of  the  conditional  mean  provided  by  the  UKF 
is  shown  to  be  correct  up  to  the  second  order  of  its  Taylor 
series  expansion.  Reference  [3]  gives  the  implementation  of 
UKF  algorithm. 

Because  the  UKF  is  the  best  in  accurately  propagating 
the  mean  and  covariance  of  the  Gaussian  approximation  to 
the  state  distribution,  it  can  be  used  to  generate  the  proposal 
distribution  for  the  particle  filter.  In  this  way,  one  obtains  a 
parametric/non-parametric  hybrid  filter  called  the 
unscented  particle  filter  (UPF). 

Ill  MAGNETIC  TARGET  TRACKING 


vectors  are  defined  in  the  same  coordinate  system,  which 
normally  are  the  sensor  coordinates.  Because  the  desired 
magnetic  moment,  ms,  is  related  to  the  ship  reference  frame, 
it  is  necessary  to  apply  a  rotation: 

S  Vy  S  Vy 

m y  -my  ,  +  mY  , 

fiTx  Vr+f  (id 

S  Vy  C  V  Y 

niy  =  m  y  ,  +  mY  . 

1  A  2  2  1  2  2 

■\IVX+VY  -\IVX+VY 

mz  -  msz 

For  a  single  sensor,  the  measurement  vector  at  time  k 
has  the  form: 

z„=(b,  B,  B;y  (12) 

These  are  the  process  and  measurement  equations  used 
by  the  Bayesian  filter  for  tracking  and  classification  of  a 
magnetic  dipole.  As  one  can  see,  the  process  function  !'(•)  in 
equation  (1)  is  linear,  and  the  measurement  function  h(»)  in 
equation  (2)  is  highly  non-linear. 


The  mathematical  model  used  for  the  target  is  a  moving 
magnetic  dipole.  The  target  is  fully  characterized  by  its 
motion  parameters  (position  and  velocity)  and  the  value  of 
the  magnetic  dipole  moment.  Its  maneuvers  and/or  a 
non-linear  trajectory  are  modeled  in  the  state  update 
through  the  process  noise  on  velocity.  The  simplifying 
assumption  made  is  that  the  target  is  moving  horizontally. 
Also  the  magnetic  mass  of  the  target  remains  constant 
during  the  measurement  and  can  be  estimated  from  the 
values  of  the  equivalent  magnetic  dipole  moment. 

Let  consider  that  the  time  increment  between  the  data 
samples  is  At  seconds,  ms  is  the  magnetic  moment  vector  of 
the  dipole  expressed  in  kA-m2,  v  is  the  velocity  vector  in 
m/sec,  and  r  is  the  position  vector  from  the  observation 
point  to  the  dipole  center  in  meters.  For  a  full 
characterization  of  the  target,  the  entire  system  at  time  step  k 
can  be  represented  by  the  state  vector: 

xk  =  (rx  rY  rz  vx  vy  mx  mY  mz  )  (8) 


The  discrete  equations  of  target  motion  are  obtained 
using  the  piece-wise  approximation: 

rx(k  )  =  rx  (k~\  )  +  Al  vx 


vx(k)=vx(k-1) 

mx  (k)  =  mx  (k  - 1) 


(9) 


and  similar  relations  exist  for  the  Y  and  Z  components  (vz  = 
0).  Tri-axial  magnetic  sensors  produce  the  observation  data. 
The  magnetic  flux  density  vector  B  at  a  given  point  due  to  a 
magnetic  dipole  is  given  by  the  formula: 
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JL 

4  n 


3(r,m)r-  r  m 
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(10) 


where  p  is  the  permeability  of  the  medium  (=  47il  0"7),  and 
<•,•>  is  the  inner  product  operator.  In  the  above  formula,  all 


TV.  SIMULATED  EXPERIMENT 


One  of  the  problems  encountered  in  designing  an 
experiment  is  the  system  observability.  A  system  is 
observable  when  its  state  vector  can  be  reconstructed  from 
the  measurements  of  its  output.  Because  the  magnetic  flux 
density  depends  on  both  states  m  and  r,  it  is  not  possible  to 
reconstruct  these  (six)  states  from  a  single  sensor,  which 
produces  at  most  3  data  points.  The  2-observers  situation 
allows  the  complete  problem  to  be  solved  [4],  i.e.  to 
determine  unambiguously  the  horizontal  position,  speed 
and  depth,  and  the  equivalent  dipole  moment  of  the  moving 
target. 


Fig.  1.  Estimated  (■)  and  true  (-)  X  and  Y  values. 

Simulated  observations  are  used  in  this  study:  a  true 
non-linear  (U-turn)  trajectory  is  assumed  for  a  magnetic 
dipole  having  a  known  moment,  observations  are  computed 
using  equation  (11)  and  (10),  and  contaminated  with 
Gaussian  noise  with  the  SNR  of  10  dB.  Two  tri-axial 
magnetometers  (observers)  placed  on  the  seafloor  have  the 


positions  given  by  Si  =  (0,  0,  0)  and  s2  =  (0,  50,  0)  meters, 
respectively.  The  target  is  represented  by  a  magnetic  dipole 
with  a  constant  moment  over  time,  ms  =  (50,  -5, 125)  kA«m2. 
The  target  approaches  the  observers  with  a  time  variable 
velocity,  and  at  a  distance  from  the  bottom  of  20m.  Discrete 
observations  (12)  are  taken  at  a  sampling  period  of  one 
second. 
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Fig.  2.  Estimated  (■)  and  true  (-)  Z  values. 


350  - 


E 

< 


O  200 

E 

o 

E 

O  150 
0 


. 

■ 

Unscented 

^article 

Filter 

/ 

N  =  36  particles 
Sensors  @(0  0 

0)  &  (0  50  0)m 

- 

60  80  100  120  140 

Time  (sec) 


Fig.  3.  Estimated  (■)  and  true  (-)  mz  values. 


In  applying  the  filtering  technique  to  the  system,  the 
initial  conditions  and  the  noise  covariance  matrixes  need  to 
be  specified.  In  the  initialization  step,  the  particles  should  be 
drawn  from  an  unknown  proposal  distribution.  The  basic 
assumption  is  that  the  target  can  approach  the  sensors  from 
any  horizontal  direction.  Therefore,  the  filter  must 
accommodate  simultaneous  alternative  hypotheses  until 
they  can  be  disambiguated  by  future  measurements.  A 
reasonable  initial  estimate  of  the  horizontal  position  is  an 
approximate  circle  around  the  sensors  with  a  radius  of  about 
200m  from  where  the  magnetic  signal  becomes  sizable.  In 
the  present  example,  36  particles  were  used  with  the 
horizontal  positions  spread  over  a  circle  every  10°  from  0° 
to  350°.  For  the  vertical  position,  an  initial  estimate  between 
zero  and  the  approximate  water  depth  can  be  given.  Because 
we  have  no  information  about  the  magnitude  of  velocity, 
acceleration  and  magnetic  dipole  moments,  a  good  initial 


estimate  of  these  vectors  are  merely  the  null  vectors. 

The  initial  covariance  matrix,  P(0|0),  gives  a  measure  of 
belief  in  the  initial  state  estimate.  It  is  assumed  that  initially 
all  the  states  are  un-correlated,  so  that  the  matrix  is  diagonal. 
This  matrix  is  not  known  and  has  to  be  sufficiently  large,  but 
the  initial  P(0|0)  is  forgotten  as  more  data  is  processed.  The 
measurement  noise  covariance  matrix  can  be  estimated 
directly  from  the  actual  data  and,  once  calculated,  it  does 
not  change  during  the  filter  run. 

The  process  noise  covariance  is  zero  for  a  deterministic 
process.  Flowever,  it  was  practically  proved  to  be  a  good 
idea  to  introduce  random  perturbations  in  the  target  position 
and  velocity.  These  small  perturbations  account  for  the 
target  maneuvers  and  prevent  divergence,  so  that  the 
process  noise  covariance  may  be  regarded  as  a  tuning 
parameter  of  the  filter. 

The  results  are  presented  in  figures  1  to  3  where  the  true 
values  of  state  variables  are  plotted  together  with  the  state 
estimates  obtained  from  this  filter.  The  performance  of  the 
UPF  is  good. 


V  CONCLUSIONS 


This  study  presents  the  application  of  the  unscented 
particle  filter  (UPF)  in  solving  the  joint  problem  of  tracking 
and  classification  of  a  target  modeled  as  an  equivalent 
magnetic  dipole  of  arbitrary  orientation.  The  problem  is 
formulated  in  state-space  form  where  the  state  variables  are 
the  position,  velocity,  and  magnetic  moment  of  the  target. 
The  advantage  of  using  the  particle  filter  for  this  application 
is  the  flexibility  in  selecting  the  initial  state  vector  to  cover 
the  possible  directions  of  arrival.  Imposing  the  initial  state 
vector  arbitrarily  represents  a  major  limitation  when  using 
the  non-linear  Kalman  filters.  On  the  other  hand,  the 
non-linear  Kalman  filter  used  in  this  study,  the  unscented 
Kalman  filter,  offers  a  better  proposal  distribution  than  the 
one  used  in  conventional  particle  filters  by  taking  into 
account  the  most  recent  measurement. 
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